Mini pruebas:
1 ciclo de asimilación dentro del experimento 20181122, para el 20181123000000 donde hay disponible observaciones del AMSU-A en NOAA15 y NOAA 18.
files <- list.files("analisis/radianzas", pattern = "amsua_n15", full.names = TRUE) %>%
.[!str_detect(., "conv")]
diag <- map(files, function(f){
meta <- unglue::unglue(f, "analisis/radianzas/asim_{sensor}_{plat}_{date}.ensmean_{exp}")
# print(f)
out <- fread(f)
# .[V10 == 1] %>%
if (file.size(f) != 0) {
out[, date := ymd_hms(meta[[1]][["date"]])] %>%
.[, exp := meta[[1]][["exp"]]]
}
out
}) %>%
rbindlist()
colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
"errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld",
"rcldp", paste0("pred", seq(8)), "date", "exp")
diag[qc == 0, .N, by = .(channel, exp)] %>%
ggplot(aes(factor(channel), exp)) +
geom_raster(aes(fill = N)) +
scale_fill_viridis_c() +
labs(x = "channel", y = "", subtitle = "QC = 0") +
theme_minimal()
diag[qc == 0 & errinv !=0, .N, by = .(channel, exp)] %>%
dcast(exp ~ channel) %>%
knitr::kable(caption = "Cantidad de observaciones potencialmente asimilables (QC == = y errinv != 0) para cada prueba y cada canal del AMSU-A")
## Using 'N' as value column. Use 'value.var' to override
| exp | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 15 |
|---|---|---|---|---|---|---|---|---|---|
| QC4new | 127 | 129 | 129 | 129 | 129 | 404 | 759 | 420 | 126 |
| QC4old | 105 | 107 | 107 | 107 | 107 | 404 | 759 | 421 | 104 |
| QC4ori | 103 | 105 | 105 | 105 | 105 | 404 | 759 | 420 | 102 |
| QC6ecmwf | 99 | 101 | 101 | 101 | 101 | 403 | 709 | 410 | 98 |
| QC6ecmwfall | 184 | 183 | 176 | 186 | 192 | 402 | 716 | 413 | 181 |
| allsky | 147 | 146 | 139 | 149 | 155 | 328 | 689 | 396 | 144 |
| clearsky | 65 | 67 | 67 | 67 | 67 | 328 | 679 | 389 | 64 |
| noQC6 | 99 | 101 | 101 | 101 | 101 | 496 | 709 | 410 | 98 |
diag[channel %in% c(1:8)] %>%
dcast(lon + lat + channel ~ exp, value.var = "qc") %>%
.[, diff := allsky - clearsky] %>%
.[diff != 0] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(color = "darkorange", size = 0.7) +
# geom_point(aes(color = errinv)) +
scale_color_viridis_c() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(.~ channel) +
labs(x = "lon",
subtitle = "Ubicación de observaciones extras en allsky respecto de clearsky") +
theme_minimal()
diag[channel %in% c(4)] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = factor(abs(qc))), size = 0.7) +
scale_color_viridis_d() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
# facet_grid(exp~channel) +
facet_wrap(~exp, ncol = 4) +
labs(subtitle = "Control de calidad para el canal 4",
x = "", y = "",
color = "QC") +
theme_minimal()
diag[channel %in% c(4:8) & qc == 0 & errinv != 0] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = tbc), size = 0.7) +
scale_color_viridis_c(option = "A") +
# scale_color_paletteer_c("ggthemes::Red-Blue Diverging") +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(channel~exp) +
labs(x = "", y = "",
color = "O-B") +
theme_minimal()
# diag[channel %in% c(1:8)] %>%
# ggplot(aes(ConvertLongitude(lon), lat)) +
# geom_point(aes(color = factor(qc))) +
# scale_color_viridis_d() +
# geom_mapa() +
# coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
# facet_grid(exp~channel) +
# theme_minimal()
# diag[channel %in% c(1:15) & errinv != 0, .N, by = .(channel, exp)] %>%
# dcast(exp~channel)
diag[, qc := abs(qc)] %>%
.[channel %in% (5), .N, by = .(exp, qc)] %>%
dcast(exp~qc) %>%
knitr::kable(caption = "CategorÃas de control de calidad por experimento para el canal 5.")
## Using 'N' as value column. Use 'value.var' to override
| exp | 0 | 3 | 4 | 8 | 50 | 51 |
|---|---|---|---|---|---|---|
| QC4new | 129 | NA | 25 | 228 | 308 | 92 |
| QC4old | 107 | NA | 23 | 138 | 307 | 207 |
| QC4ori | 105 | NA | 21 | 139 | 308 | 209 |
| QC6ecmwf | 101 | NA | 18 | 142 | 310 | 211 |
| QC6ecmwfall | 192 | 3 | 21 | 87 | 311 | 168 |
| allsky | 155 | 3 | 4 | 38 | 450 | 132 |
| clearsky | 67 | NA | 3 | 91 | 450 | 171 |
| noQC6 | 101 | NA | 34 | 190 | NA | 457 |
diag[channel %in% c(1:15) & qc == 0 & errinv > 0.000031623,
.N, by = .(exp)] %>%
dcast(exp~.) %>%
.[order(-.)] %>%
knitr::kable(caption = "Cantidad total de observaciones potencialmente asimilables")
## Using 'N' as value column. Use 'value.var' to override
| exp | . |
|---|---|
| QC6ecmwfall | 2633 |
| QC4new | 2352 |
| allsky | 2293 |
| QC4old | 2221 |
| noQC6 | 2216 |
| QC4ori | 2208 |
| QC6ecmwf | 2123 |
| clearsky | 1793 |
diag[channel %in% c(4:8)] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = errinv)) +
scale_color_viridis_c() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(channel~exp) +
labs(x = "", y = "") +
theme_minimal()
files <- list.files("analisis/radianzas", pattern = "analysis", full.names = TRUE) %>%
.[!str_detect(., "asim")]
update <- map(files, function(f) {
descriptores <- unglue::unglue(f, c("analisis/radianzas/analysis.ensmean_{exp}"))
ana <- ReadNetCDF(f, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR",
lon = "XLONG", lat = "XLAT")) %>%
.[, p := p + PB] %>%
.[, t := tk(t, p)] %>%
.[, rh := rh(qv, p, t)] %>%
.[, td := td(qv, p) + 273.15] %>%
.[, ":="(Time = NULL,
# west_east = NULL,
# south_north = NULL,
qv = NULL,
PB = NULL)] %>%
.[, c("u", "v") := uvmet(f)] %>%
.[, ":="(exp = descriptores[[1]][["exp"]])] %>%
.[]
}) %>%
rbindlist()
guess <- ReadNetCDF("analisis/radianzas/wrfarw.ensmean",
vars = c(p = "P", "PB", t = "T", qv = "QVAPOR",
lon = "XLONG", lat = "XLAT")) %>%
.[, p := p + PB] %>%
.[, t := tk(t, p)] %>%
.[, rh := rh(qv, p, t)] %>%
.[, td := td(qv, p) + 273.15] %>%
.[, ":="(Time = NULL,
# west_east = NULL,
# south_north = NULL,
qv = NULL,
PB = NULL)] %>%
.[, c("u", "v") := uvmet("analisis/radianzas/wrfarw.ensmean")] %>%
# .[, ":="(exp = "guess")] %>%
.[]
update <- guess[update, on = c("lon", "lat", "bottom_top", "south_north", "west_east")]
update[, .(mean_ana = mean(i.t, na.rm = TRUE),
mean_guess = mean(t, na.rm = TRUE)), by = .(bottom_top, west_east, exp)] %>%
.[exp %in% c("QC4new", "QC4old")] %>%
ggplot(aes(west_east, bottom_top)) +
geom_contour_fill(aes(z = (mean_ana - mean_guess))) +
scale_fill_divergent() +
labs(fill = "update T") +
facet_wrap(~ exp) +
theme_minimal()
update[, .(mean_ana = mean(i.t, na.rm = TRUE),
mean_guess = mean(t, na.rm = TRUE)), by = .(bottom_top, south_north, exp)] %>%
.[exp %in% c("QC4new", "QC4old")] %>%
ggplot(aes(south_north, bottom_top)) +
geom_contour_fill(aes(z = (mean_ana - mean_guess))) +
scale_fill_divergent() +
labs(fill = "update T") +
facet_wrap(~ exp) +
theme_minimal()
La diferencia entre ambas pruebas para la temperatura es en promedio del ordende 10-3 o 10-4 y un poco menor en los niveles mas altos. El rango de la diferencia está aproximadamente entre -0.4 y 0.2 K.
update %>%
dcast(bottom_top + south_north + west_east ~ exp, value.var = "i.t") %>%
.[, .(diff = mean(abs(QC4new - QC4old)),
max = max(QC4new - QC4old),
min = min(QC4new - QC4old)), by = .(bottom_top)] %>%
melt(measure.vars = c("diff", "max", "min")) %>%
ggplot(aes(value, bottom_top)) +
geom_vline(xintercept = 0, color = "darkgray") +
geom_path(aes(color = variable)) +
scale_color_viridis_d() +
labs(y = "sigma level", color = "") +
theme_minimal()
files <- list.files("analisis/diagfiles/E7", pattern = "asim", full.names = TRUE) %>%
.[!str_detect(., "conv")]
diag <- map(files, function(f){
meta <- unglue::unglue(f, "analisis/diagfiles/{exp}/asim_{sensor}_{plat}_{date}.ensmean")
# print(f)
out <- fread(f)
# .[V10 == 1] %>%
if (file.size(f) != 0) {
out[, date := ymd_hms(meta[[1]][["date"]])]
}
out
}) %>%
rbindlist()
colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
"errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld",
"rcldp", paste0("pred", seq(8)), "date")
satinfo <- fread("global_satinfo.txt") %>%
setnames(old = c("!sensor/instr/sat", "chan"), new = c("sensor", "channel"))
diag[, .N, by = .(sensor, channel, qc)] %>%
satinfo[., on = c("sensor", "channel")] %>%
ggplot(aes(factor(channel), factor(abs(qc)))) +
geom_tile(aes(fill = N)) +
scale_fill_viridis_c() +
geom_point(data = ~.x[iuse != 1], shape = 4, color = "grey10") +
# scale_shape_manual(values = c(0, 4)) +
facet_wrap(~sensor, scales = "free_x") +
labs(x = "channel", y = "QC",
caption = "Con 'x' se marcan los canales rechazados por satinfo") +
theme_minimal()